

## ===================== Input tcl File ===================================
inputfiles = ["/home/chenting/Chenting/bin/FrameTrial/NodeEles1.tcl", "/home/chenting/Chenting/bin/FrameTrial/NodeEles2.tcl", "/home/chenting/Chenting/bin/FrameTrial/BeamElements.tcl"]

### ====================== Do you wish to display result? ==================
DispRes = True
StressRes = False
PlotOriginal = False
nframe = 299
xscale = 5
yscale = 5
zscale = 5

if DispRes:
    resfile = "/home/chenting/Chenting/bin/FrameTrial/NodeDisp.out"

if StressRes:
    ShellNLResfile = "/home/chenting/Chenting/bin/DualPlateRCorner/EleStress.out"
    E = 219000000000
    pos = 0.3
    fy = 332000000
    z = 0.003
    h = 0.006
    IPFactor = 0.774597

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.gca(projection='3d')

NodeList=[]
FrameEleList = []
ShellEleList = []

def pprint(A):
    n = len(A)
    for i in range(0, n):
        line = ""
        for j in range(0, n+1):
            line += str(A[i][j]) + "\t"
            if j == n-1:
                line += "| "
        print(line)
    print("")


def myGauss(A):
    n = len(A)

    for i in range(0, n):
        # Search for maximum in this column
        maxEl = abs(float(A[i][i]))
        maxRow = i
        for k in range(i+1, n):
            if abs(A[k][i]) > maxEl:
                maxEl = abs(float(A[k][i]))
                maxRow = k

        # Swap maximum row with current row (column by column)
        for k in range(i, n+1):
            tmp = float(A[maxRow][k])
            A[maxRow][k] = float(A[i][k])
            A[i][k] = tmp

        # Make all rows below this one 0 in current column
        for k in range(i+1, n):
            c = -float(A[k][i])/float(A[i][i])
            for j in range(i, n+1):
                if i == j:
                    A[k][j] = 0
                else:
                    A[k][j] += c * float(A[i][j])

    # Solve equation Ax=b for an upper triangular matrix A
    x = [0 for i in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = float(A[i][n])/A[i][i]
        for k in range(i-1, -1, -1):
            A[k][n] -= A[k][i] * x[i]
    return x

class Node:
    def __init__(self, n, x, y, z):
        self.n = n
        self.x = float(x)
        self.y = float(y)
        self.z = float(z)

    def x(self):
        return self.x

    def y(self):
        return self.y

    def z(self):
        return self.z

    def n(self):
        return self.n


class FrameEle:
    def __init__(self, type, n, n1, n2):
        self.type = type
        self.n = n
        self.n1 = n1
        self.n2 = n2

    def type(self):
        return self.type
    def n(self):
        return self.n
    def n1(self):
        return self.n1
    def n2(self):
        return self.n2



class ShellEle:
    def __init__(self, type, n, n1, n2, n3, n4, n5=0, n6=0, n7=0, n8=0, n9=0):
        self.type = type
        self.n = n
        self.n1 = n1
        self.n2 = n2
        self.n3 = n3
        self.n4 = n4
        self.n5 = n5
        self.n6 = n6
        self.n7 = n7
        self.n8 = n8
        self.n9 = n9

    def type(self):
        return self.type
    def n(self):
        return self.n
    def n1(self):
        return self.n1
    def n2(self):
        return self.n2
    def n3(self):
        return self.n3
    def n4(self):
        return self.n4
    def n5(self):
        return self.n5
    def n6(self):
        return self.n6
    def n7(self):
        return self.n7
    def n8(self):
        return self.n8
    def n9(self):
        return self.n9

for iter in range(0, len(inputfiles)):
    content=open(inputfiles[iter])
    for line in content:
        splited = line.split()
        if len(splited) > 0:
            if splited[0] == 'node':
                NodeList.append(Node(splited[1], splited[2], splited[3], splited[4]))

            if splited[0] == 'element':
                if splited[1] == 'ShellMITC4' or splited[1] == 'ShellNLDKGQ' or splited[1] == 'ShellDKGQ':
                    ShellEleList.append(ShellEle(splited[1], splited[2], splited[3], splited[4], splited[5], splited[6]))
                if splited[1] == 'ShellNL':
                    ShellEleList.append(ShellEle(splited[1], splited[2], splited[3], splited[4], splited[5], splited[6], splited[7], splited[8], splited[9], splited[10], splited[11]))
                if splited[1] == 'forceBeamColumn' or splited[1] == 'dispBeamColumn' or splited[1] == 'elasticBeamColumn' or splited[1] == 'beamWithHinges':
                    FrameEleList.append(FrameEle(splited[1], splited[2], splited[3], splited[4]))
    content.close()

nNodes = len(NodeList)
nShells = len(ShellEleList)
nBeams = len(FrameEleList)
print('a total of %d nodes'%nNodes)

NodeList2 = sorted(NodeList, key = lambda i: int(i.n))
ShellEleList2 = sorted(ShellEleList, key = lambda i: int(i.n))
FrameEleList2 = sorted(FrameEleList, key = lambda i: int(i.n))

xs = map(lambda i: NodeList2[i].x, range(0, nNodes))
ys = map(lambda i: NodeList2[i].y, range(0, nNodes))
zs = map(lambda i: NodeList2[i].z, range(0, nNodes))

## ====== Plot Nodes ==============================
if PlotOriginal:
    ax.scatter(xs, ys, zs)

for i in range(0, nShells):
    NodeI = filter(lambda j: NodeList2[j].n == ShellEleList[i].n1, range(0, nNodes))[0]
    NodeJ = filter(lambda j: NodeList2[j].n == ShellEleList[i].n2, range(0, nNodes))[0]
    NodeK = filter(lambda j: NodeList2[j].n == ShellEleList[i].n3, range(0, nNodes))[0]
    NodeL = filter(lambda j: NodeList2[j].n == ShellEleList[i].n4, range(0, nNodes))[0]
    NodeQ = filter(lambda j: NodeList2[j].n == ShellEleList[i].n9, range(0, nNodes))[0]
    x1 = NodeList2[NodeI].x
    x2 = NodeList2[NodeJ].x
    x3 = NodeList2[NodeK].x
    x4 = NodeList2[NodeL].x
    y1 = NodeList2[NodeI].y
    y2 = NodeList2[NodeJ].y
    y3 = NodeList2[NodeK].y
    y4 = NodeList2[NodeL].y
    z1 = NodeList2[NodeI].z
    z2 = NodeList2[NodeJ].z
    z3 = NodeList2[NodeK].z
    z4 = NodeList2[NodeL].z
    x9 = NodeList2[NodeQ].x
    y9 = NodeList2[NodeQ].y
    z9 = NodeList2[NodeQ].z
    if PlotOriginal:
        for j in range(1, 4):
            eval("ax.plot([x%d, x%d], [y%d, y%d], [z%d, z%d], 'r-')"%(j, j+1, j, j+1, j, j+1))
        ax.plot([x4, x1], [y4, y1], [z4, z1], 'r-')
        ax.plot([x9, x1], [y9, y1], [z9, z1], 'r-')


nline = 0
xlist = filter(lambda i: i%6 == 1, range(0, nNodes*6+1))
ylist = filter(lambda i: i%6 == 2, range(0, nNodes*6+1))
zlist = filter(lambda i: i%6 == 3, range(0, nNodes*6+1))
dx = []
dy = []
dz = []

if DispRes:
    content = open(resfile)
    for line in content:
        nline = nline + 1
        if nline == nframe:
            splited = line.split()
            for obj in xlist:
                dx.append(float(splited[obj]))
            for obj in ylist:
                dy.append(float(splited[obj]))
            for obj in zlist:
                dz.append(float(splited[obj]))
            break

    print(len(xs), len(dx))
    print('a total of %d nodes for resualts'%len(dx))
    resx = map(lambda i: xs[i] + xscale * dx[i], range(0, nNodes))
    resy = map(lambda i: ys[i] + yscale * dy[i], range(0, nNodes))
    resz = map(lambda i: zs[i] + zscale * dz[i], range(0, nNodes))

    ax.scatter(resx, resy, resz)

    for i in range(0, nShells):
        if ShellEleList2[i].type == 'ShellNL':
            NodeI = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n1, range(0, nNodes))[0]
            NodeJ = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n2, range(0, nNodes))[0]
            NodeK = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n3, range(0, nNodes))[0]
            NodeL = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n4, range(0, nNodes))[0]
            NodeM = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n5, range(0, nNodes))[0]
            NodeN = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n6, range(0, nNodes))[0]
            NodeO = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n7, range(0, nNodes))[0]
            NodeP = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n8, range(0, nNodes))[0]
            NodeQ = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n9, range(0, nNodes))[0]
            x1 = resx[NodeI]
            x2 = resx[NodeJ]
            x3 = resx[NodeK]
            x4 = resx[NodeL]
            x5 = resx[NodeM]
            x6 = resx[NodeN]
            x7 = resx[NodeO]
            x8 = resx[NodeP]
            x9 = resx[NodeQ]

            y1 = resy[NodeI]
            y2 = resy[NodeJ]
            y3 = resy[NodeK]
            y4 = resy[NodeL]
            y5 = resy[NodeM]
            y6 = resy[NodeN]
            y7 = resy[NodeO]
            y8 = resy[NodeP]
            y9 = resy[NodeQ]

            z1 = resz[NodeI]
            z2 = resz[NodeJ]
            z3 = resz[NodeK]
            z4 = resz[NodeL]
            z5 = resz[NodeM]
            z6 = resz[NodeN]
            z7 = resz[NodeO]
            z8 = resz[NodeP]
            z9 = resz[NodeQ]
            for j in range(1, 4):
                eval("ax.plot([x%d, x%d], [y%d, y%d], [z%d, z%d], 'b-')" % (j, j + 1, j, j + 1, j, j + 1))
            ax.plot([x4, x1], [y4, y1], [z4, z1], 'b-')

        if ShellEleList2[i].type == 'ShellMITC4' or ShellEleList2[i].type == 'ShellNLDKGQ':
            NodeI = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n1, range(0, nNodes))[0]
            NodeJ = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n2, range(0, nNodes))[0]
            NodeK = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n3, range(0, nNodes))[0]
            NodeL = filter(lambda j: NodeList2[j].n == ShellEleList2[i].n4, range(0, nNodes))[0]

            x1 = resx[NodeI]
            x2 = resx[NodeJ]
            x3 = resx[NodeK]
            x4 = resx[NodeL]

            y1 = resy[NodeI]
            y2 = resy[NodeJ]
            y3 = resy[NodeK]
            y4 = resy[NodeL]

            z1 = resz[NodeI]
            z2 = resz[NodeJ]
            z3 = resz[NodeK]
            z4 = resz[NodeL]

            for j in range(1, 4):
                eval("ax.plot([x%d, x%d], [y%d, y%d], [z%d, z%d], 'b-')" % (j, j + 1, j, j + 1, j, j + 1))
            ax.plot([x4, x1], [y4, y1], [z4, z1], 'b-')

    for i in range(0, nBeams):
        if FrameEleList2[i].type == 'forceBeamColumn' or FrameEleList2[i].type == 'dispBeamColumn' or FrameEleList2[i].type == 'elasticBeamColumn' or FrameEleList2[i].type == 'beamWithHinges':
            NodeI = filter(lambda j: NodeList2[j].n == FrameEleList2[i].n1, range(0, nNodes))[0]
            NodeJ = filter(lambda j: NodeList2[j].n == FrameEleList2[i].n2, range(0, nNodes))[0]

            x1 = resx[NodeI]
            x2 = resx[NodeJ]

            y1 = resy[NodeI]
            y2 = resy[NodeJ]

            z1 = resz[NodeI]
            z2 = resz[NodeJ]

            ax.plot([x1, x2], [y1, y2], [z1, z2], 'b-')

    content.close()

if StressRes:
    for i in range(1, 9):
        exec('s%dlist = []' % i)
    for i in range(1, 10):
        s1list.append(filter(lambda j: j % 84 == 1 + 8*(i - 1), range(0, nShells * 84 + 1)))
        s2list.append(filter(lambda j: j % 84 == 2 + 8*(i - 1), range(0, nShells * 84 + 1)))
        s3list.append(filter(lambda j: j % 84 == 3 + 8*(i - 1), range(0, nShells * 84 + 1)))
        s4list.append(filter(lambda j: j % 84 == 4 + 8*(i - 1), range(0, nShells * 84 + 1)))
        s5list.append(filter(lambda j: j % 84 == 5 + 8*(i - 1), range(0, nShells * 84 + 1)))
        s6list.append(filter(lambda j: j % 84 == 6 + 8*(i - 1), range(0, nShells * 84 + 1)))
        s7list.append(filter(lambda j: j % 84 == 7 + 8*(i - 1), range(0, nShells * 84 + 1)))
        s8list.append(filter(lambda j: j % 84 == 8 + 8*(i - 1), range(0, nShells * 84 + 1)))

    content = open(ShellNLResfile)
    for i in range(1, 9):
        exec('s%dv = [[],[],[],[],[],[],[],[],[]]' % i)
    for line in content:
        nline = nline + 1
        if nline == nframe:
            splited = line.split()
            print(len(splited))
            for j in range(0, 9):
                for obj in s1list[j]:
                    s1v[j].append(float(splited[obj]))
                for obj in s2list[j]:
                    s2v[j].append(float(splited[obj]))
                for obj in s3list[j]:
                    s3v[j].append(float(splited[obj]))
                for obj in s4list[j]:
                    s4v[j].append(float(splited[obj]))
                for obj in s5list[j]:
                    s5v[j].append(float(splited[obj]))
                for obj in s6list[j]:
                    s6v[j].append(float(splited[obj]))
                for obj in s7list[j]:
                    s7v[j].append(float(splited[obj]))
                for obj in s8list[j]:
                    s8v[j].append(float(splited[obj]))

    eletrans = []

    for i in range(10000, 20000):
        temp = filter(lambda j: int(ShellEleList[j].n) == i, range(0, nShells))
        if len(temp) == 1:
            eletrans.append(temp[0])
    k = 0

    VMStressMax = 0

    for i in eletrans:
        if ShellEleList[i].type == 'ShellNL':
            NodeI = filter(lambda j: NodeList2[j].n == ShellEleList[i].n1, range(0, nNodes))[0]
            NodeJ = filter(lambda j: NodeList2[j].n == ShellEleList[i].n2, range(0, nNodes))[0]
            NodeK = filter(lambda j: NodeList2[j].n == ShellEleList[i].n3, range(0, nNodes))[0]
            NodeL = filter(lambda j: NodeList2[j].n == ShellEleList[i].n4, range(0, nNodes))[0]
            NodeM = filter(lambda j: NodeList2[j].n == ShellEleList[i].n5, range(0, nNodes))[0]
            NodeN = filter(lambda j: NodeList2[j].n == ShellEleList[i].n6, range(0, nNodes))[0]
            NodeO = filter(lambda j: NodeList2[j].n == ShellEleList[i].n7, range(0, nNodes))[0]
            NodeP = filter(lambda j: NodeList2[j].n == ShellEleList[i].n8, range(0, nNodes))[0]
            NodeQ = filter(lambda j: NodeList2[j].n == ShellEleList[i].n9, range(0, nNodes))[0]
            x1 = NodeList2[NodeI].x
            x2 = NodeList2[NodeJ].x
            x3 = NodeList2[NodeK].x
            x4 = NodeList2[NodeL].x
            y1 = NodeList2[NodeI].y
            y2 = NodeList2[NodeJ].y
            y3 = NodeList2[NodeK].y
            y4 = NodeList2[NodeL].y
            z1 = NodeList2[NodeI].z
            z2 = NodeList2[NodeJ].z
            z3 = NodeList2[NodeK].z
            z4 = NodeList2[NodeL].z

            #### ============================== Need to find which direction is perpendiculer to the shell surface ======================
            alplist = ['I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q']

            for j in range(0, 9):
                for iter in range(1, 9):
                    exec('s%dv[j][k] = (s%dv[j][k] - s%dv[8][k])/IPFactor + s%dv[8][k]' % (iter, iter, iter, iter))
                #print('The %d th element, the %d th node, %f, %f, %f, %f, %f, %f, %f, %f' % (i, j, s1v[j][k], s2v[j][k], s3v[j][k], s4v[j][k], s5v[j][k], s6v[j][k], s7v[j][k], s8v[j][k]))
                b = s1v[j][k]/h
                a = 12*s4v[j][k]/(h**3)
                stressx = a*z + b

                b = s2v[j][k] / h
                a = 12 * s5v[j][k] / (h ** 3)
                stressy = a * z + b

                b = s3v[j][k] / h
                a = 12 * s6v[j][k] / (h ** 3)
                stressxy = a * z + b

                stressxz = s7v[j][k] / h
                stressyz = s8v[j][k] / h

                VMStress = (stressx - stressy)**2 + stressy**2 + stressx**2 + 6*(stressxy**2 + stressxz**2 + stressyz**2)
                VMStress = (VMStress*0.5)**0.5

                if VMStress > VMStressMax:
                    VMStressMax = VMStress

                print(ShellEleList[i].n)

                if int(ShellEleList[i].n) > 20000:
                    VMStress = -VMStress

                print('The %d th element, the %d th node, %f' % (i+1, j+1, VMStress))
                IPOrderlist = [1, 5, 2, 6, 3, 7, 4, 8, 9]
                exec('VMStress%d = VMStress' % IPOrderlist[j])

            k = k + 1

            for iter in range(1, 4):
                zscale = 1
                exec("ax.plot([x%d, x%d], [y%d, y%d], [z%d + VMStress%d/zscale, z%d + VMStress%d/zscale], 'b-')" % (iter, iter + 1, iter, iter + 1, iter, iter, iter + 1, iter + 1))
                ax.plot([x4, x1], [y4, y1], [z4 + VMStress4/zscale, z1 + VMStress1/zscale], 'b-')

    print('Max Von Mise Stress is %f' % VMStressMax)
    content.close()

ax.auto_scale_xyz([0, 6], [-3, 3], [0, 6])
#ax.auto_scale_xyz([-0.2, 0.2], [-0.2, 0.2], [5.8, 6.2])
plt.show()

